rm(list = ls())
options(warn=-1)
library(rgl)
library(misc3d)
knitr::knit_hooks$set(webgl = hook_webgl)

#molecular-atlas-master is available at https://github.com/cantin-ortiz/molecular-atlas
path.bin <- './molecular-atlas-master/bin'
path.matrices <-'./molecular-atlas-master/data'

load(paste(path.matrices ,'atlasspots.RData',sep='/'))
load(paste(path.matrices , 'vivid-colors.RData', sep='/'))

if(!exists('aligned.atlas'))
  load(paste(path.matrices ,'alignedAtlas.RData',sep='/'))

if(!exists('atlas.stereo'))
  load(paste(path.matrices ,'atlasStereo.RData',sep='/'))

source(paste(path.bin,'plotFunctions.R',sep='/'))
source(paste(path.bin,'execFunctions.R',sep='/'))
source(paste(path.bin,'araAtlasFunctions.R',sep='/'))
source(paste(path.bin,'smoothingFunctions.R',sep='/'))
source(paste(path.bin,'layerFunctions.R',sep='/'))
source(paste(path.bin,'icFunctions.R',sep='/'))
source(paste(path.bin,'tsneFunctions.R',sep='/'))
source(paste(path.bin,'similarityIndexFunctions.R',sep='/'))
source(paste(path.bin,'scRemapFunctions.R',sep='/'))
source(paste(path.bin,'constantParameters.R',sep='/'))
source(paste(path.bin,'allenAnnotationsFunctions.R',sep='/'))
source(paste(path.bin,'meshCutsFunctions.R',sep='/'))

Load results and 3D locations of spots.

#STitch3D's results and CCFv3 coordinates.
load("spotstable.RData")
load("VOLUMESMALL.RData")
load("VOLUME.RData")
for (i in (0:34)){
  prop <- read.table(paste0("./res_mouse_brain/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
library(stringr)
prop_all$X <- str_split_fixed(prop_all$X, "-", 2)[, 1]
spots.table <- spots.table[prop_all$X, ]
spots.table$X <- rownames(spots.table)
spots.table <- merge(spots.table, prop_all, by=c("X"))
spots.table.model <- spots.table

#STitch3D's results and constructed 3D locations of spots.
load("spotstable.RData")
for (i in (0:34)){
  prop <- read.table(paste0("./res_mouse_brain/prop_slice", i, ".csv"), sep=",", header=TRUE)
  if (i == 0){
    prop_all <- prop
  }else{
    prop_all <- rbind(prop_all, prop)
  }
}
prop_all$X <- str_split_fixed(prop_all$X, "-", 2)[, 1]
spots.table <- spots.table[prop_all$X, ]
spots.table$X <- rownames(spots.table)
spots.table <- merge(spots.table, prop_all, by=c("X"))
colnames(spots.table)[10] <- "x_pixel"
colnames(spots.table)[11] <- "y_pixel"
coor_3d <- read.table(paste0("./res_mouse_brain/3D_coordinates.csv"), sep=",", header=TRUE)
coor_3d$X <- str_split_fixed(coor_3d$X, "-", 2)[, 1]
spots.table <- merge(spots.table, coor_3d, by=c("X"))
spots.table.ours <- spots.table

#clustering results
clusters <- read.table(paste0("./res_mouse_brain/clustering_result.csv"), sep=",", header=TRUE)
clusters$X <- str_split_fixed(clusters$X, "-", 2)[, 1]

louvain_colors <- c('#8c564b', '#ff7f0e', '#2ca02c', '#17becf', 
                    '#9467bd', '#1f77b4', '#e377c2', '#7f7f7f', 
                    '#bcbd22', '#d62728', '#aec7e8')
hpc_colors = c("#d62728", "#FFFF00", "#1f77b4", "#2ca02c")

STitch3D’s spatial domain detection result reveals layer structures in cortex.

spots.table <- spots.table.ours
spots.table <- merge(spots.table, clusters, by=c("X"))

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   1

par3d(persp)

## NULL

um <- c(-0.03762093, -0.06319845, 0.997291744, 0,
        0.84057891, -0.54168296, -0.002617532, 0,
        0.54038119, 0.83820403, 0.073501527, 0,
        0, 0, 0, 1)
view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
  
spheres3d(spots.table[spots.table$louvain==1, ]$x, 
          spots.table[spots.table$louvain==1, ]$y, 
          spots.table[spots.table$louvain==1, ]$z,
          col = louvain_colors[2], radius=0.5, 
          alpha=1)
spheres3d(spots.table[spots.table$louvain==2, ]$x, 
          spots.table[spots.table$louvain==2, ]$y, 
          spots.table[spots.table$louvain==2, ]$z,
          col = louvain_colors[3], radius=0.5, 
          alpha=1)
spheres3d(spots.table[spots.table$louvain==5, ]$x, 
          spots.table[spots.table$louvain==5, ]$y, 
          spots.table[spots.table$louvain==5, ]$z,
          col = louvain_colors[6], radius=0.5, 
          alpha=1)
spheres3d(spots.table$x, spots.table$y, spots.table$z,
          col = 'gray', radius=0.5, 
          alpha=0.02)
decorate3d()

spots.table <- spots.table.model
spots.table <- merge(spots.table, clusters, by=c("X"))

open3d(windowRect = c(0, 0, 720, 720))

## glX 
##   3

par3d(persp)

## NULL

um <- c(0.95741981, -0.2865749, -0.03495293, 0,
        -0.07430246, -0.1276060, -0.98903739, 0,
        0.27897382, 0.9495215, -0.14346580, 0,
        0, 0, 0, 1)
view3d(userMatrix = matrix(um, byrow=TRUE, nrow=4))
drawScene.rgl(list(VOLUMESMALL))
smp.AP <- 0
smp.ML <- 0
spheres3d(spots.table[spots.table$louvain==1, ]$AP.paxTOallen - 530/2 + smp.AP,
          -spots.table[spots.table$louvain==1, ]$DV * 1000/25 - 320/2,
          spots.table[spots.table$louvain==1, ]$ML * 1000/25 + smp.ML,
          col = louvain_colors[2], radius=5, alpha=1)
spheres3d(spots.table[spots.table$louvain==1, ]$AP.paxTOallen - 530/2 + smp.AP,
          -spots.table[spots.table$louvain==2, ]$DV * 1000/25 - 320/2,
          spots.table[spots.table$louvain==2, ]$ML * 1000/25 + smp.ML,
          col = louvain_colors[3], radius=5, alpha=1)
spheres3d(spots.table[spots.table$louvain==1, ]$AP.paxTOallen - 530/2 + smp.AP,
          -spots.table[spots.table$louvain==5, ]$DV * 1000/25 - 320/2,
          spots.table[spots.table$louvain==5, ]$ML * 1000/25 + smp.ML,
          col = louvain_colors[6], radius=5, alpha=1)